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The postsynaptic potentials of pyramidal neurons have a non-Gaussian amplitude 
distribution with a heavy tail in both hippocampus and neocortex. Such distributions 
of synaptic weights were recently shown to generate spontaneous internal noise 
optimal for spike propagation in recurrent cortical circuits. However, whether this internal 
noise generation by heavy-tailed weight distributions is possible for and beneficial to 
other computational functions remains unknown. To clarify this point, we construct an 
associative memory (AM) network model of spiking neurons that stores multiple memory 
patterns in a connection matrix with a lognormal weight distribution. In AM networks, 
non-retrieved memory patterns generate a cross-talk noise that severely disturbs memory 
recall. We demonstrate that neurons encoding a retrieved memory pattern and those 
encoding non-retrieved memory patterns have different subthreshold membrane-potential 
distributions in our model. Consequently, the probability of responding to inputs at strong 
synapses increases for the encoding neurons, whereas it decreases for the non-encoding 
neurons. Our results imply that heavy-tailed distributions of connection weights can 
generate noise useful for AM recall. 
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INTRODUCTION 

The organization of neuronal wiring determines the flow of 
information in neural circuits and hence has significant impli- 
cations for functions of the circuits. A number of recent studies 
revealed the non-random features of neuronal wiring in cortical 
circuits (Markram, 1997; Kalisman et al., 2005; Song et al., 2005; 
Yoshimura et al., 2005; Koulakov et al, 2009; Lefort et al, 2009; 
Yassin et al., 2010; Perin et al., 2011). These features include not 
only the complex topology of neuronal wiring, but also the non- 
Gaussian nature of the distributions of amplitudes of excitatory 
postsynaptic potentials (EPSPs), which is typically represented by 
the presence of long tails (also called "heavy tails") in the distri- 
butions. In fact, the amplitude distributions of EPSPs are well 
described by lognormal distributions in both neocortex (Song 
et al., 2005; Lefort et al., 2009) and hippocampus (Ikegaya et al., 
2013). In reality, the amplitude of EPSP between cortical neu- 
rons may represent the total strength of multiple synaptic contacts 
made by one of the neurons on the other. However, hereafter we 
simply call the EPSP amplitude the "synaptic weight" from one 
neuron to the other. 

A lognormal weight distribution implies that a small number 
of very strong connections are present in local cortical circuits 
and carry a large amount of the total weight on a cortical neu- 
ron, while the majority of synapses are weak (Holmgren et al., 
2003; Binzegger et al., 2004). Such EPSP-amplitude distributions 
of AMPA receptor-mediated synapses significantly influence the 



dynamical properties of stable network states (Koulakov et al., 
2009; Roxin et al., 2011). In particular, we have recently shown 
that spontaneous cortical activity emerges from the coopera- 
tion between strong-sparse and weak-dense (SSWD) synapses in 
a heavy-tailed EPSP-amplitude distribution of AMPA recurrent 
synapses (Teramae et al., 2012). Such a recurrent network can 
generate internal noise optimal for stochastic resonance effects 
on spike-based communications between neurons. Moreover, 
lognormally-connected recurrent networks combined with highly 
non-random properties of cortical synaptic connections (Prill 
et al., 2005; Song et al., 2005; Perin et al., 2011) generate working 
memory-like bistable network states (Klinshov et al, unpub- 
lished data). These results seem to indicate that heavy- tailed 
weight distributions play active roles in stochastic cortical infor- 
mation processing including the pattern-recall operation in the 
hippocampus. However, these roles remain largely unknown. 

Hippocampal CA3 has sparse recurrent excitatory connec- 
tions and is thought to perform a pattern completion oper- 
ation in memory recall (Nakazawa et al., 2002). Properties 
of pattern completion, such as storage capacity and memory 
retrieval dynamics, have been extensively studied in the statis- 
tical mechanics of Hopfield associative memory (AM) network 
and its variants of formal neurons with binary or analog out- 
puts (Hopfield, 1982, 1984; Amit et al., 1985; Derrida et al, 1987; 
Tsodyks and Feigel'man, 1988; Shiino and Fukai, 1992; Coolen 
and Sherrington, 1993; Okada, 1995). AM models of spiking 
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neurons have also been studied although they generally exhibit 
a much poorer ability for memory storage than models of binary 
or analog neurons (Gerstner and Hemmen, 1992; Lansner and 
Fransen, 1992; Treves, 1993; Amit and Brunei, 1997; Maass and 
Natschlager, 1997; Sommer and Wennekers, 2001; Curti et al., 
2004; Latham and Nirenberg, 2004; Aviel et al, 2005; Roudi and 
Latham, 2007). In most of the previous models, the weights of 
synaptic connections typically obey a Gaussian distribution, par- 
ticularly when the number of embedded patterns is extensively 
large. However, since the amplitudes of EPSPs between pyrami- 
dal cells were shown to obey a heavy-tailed distribution in CA3 
(Ikegaya et al., 2013), here we construct AM network models of 
spiking neurons having lognormal weight distributions of recur- 
rent connections and investigate possible implications of such 
non-Gaussian distributions in the memory retrieval dynamics. 
In particular, we explore the possibility that such a network can 
generate internal noise useful for the retrieval of an embedded 
memory pattern. 

MATERIALS AND METHODS 

We developed a recurrent network model of integrate-and-fire 
neurons that stores information on a multiple number of memory 
patterns in the weights of excitatory synaptic connections obeying 
a lognormal distribution. Extremely strong synapses are rare in 
such a weight distribution and the majority of synapses are weak. 
Therefore, we may term neural networks with long- tailed synaptic 
weight distributions "SSWD networks." 

NETWORK DYNAMICS 

The model consists of Nb excitatory neurons and Nj inhibitory 
neurons. The excitatory neurons are connected with each other by 
a lognormally weight-distributed Hebbian-type synaptic matrix 



that is generated from random memory patterns, and each 
inhibitory neuron is connected randomly with the excitatory neu- 
rons and the other inhibitory neurons (Figure 1 A). The neural 
dynamics in this recurrent network is described by conductance- 
based leaky integrate-and-fire neurons as: 

^ = - — (Vi - Vl) - gf (v, - V E ) - g\ (v, - V;) (1) 

(XX "^m 

where v; is the membrane potential of neuron i, x m is the 
membrane time constant, and Ve> Vi, and Vi are the rever- 
sal potentials of AMPA-receptor-mediated excitatory synaptic 
current, inhibitory synaptic current, and leak current, respec- 
tively. The conductances of excitatory and inhibitory synapses are 
measured in units of [time -1 ] and obey. 

i s i 

(2) 

where xf is the decay constant of excitatory or inhibitory synap- 
tic current, s- ' is the spike times of excitatory or inhibitory 
neuron j, df Y is delay from neuron j to i, and the element of 
the adjacent matrix c^ Y is 1 if neurons j and i are connected 
and otherwise 0. Indices X and Y indicate whether pre- or post- 
neurons are excitatory or inhibitory. The maximum conductance 
G EE of excitatory-to-excitatory synapses was determined such 
that the amplitudes of their EPSPs obey a lognormal distribu- 
tion lnN(|i, cr 2 ), with parameters u, and a 2 being the mean and 
variance of the normal distribution. Other types of conductance 
(E-to-I, I-to-E, J-to-7) were fixed at constant values. 




Vij[mV] 




FIGURE 1 | The construction of an associative memory network with 
strong-sparse and weak-dense synapses. (A) Schematic illustration of 
strong-sparse weak-dense connected network. (B) The cumulative 
distribution function (cdf) was calculated from the probability distribution of 
the elements of the Hebbian connection matrix Jn [left). Then, the EPSP 



amplitude 14 of this synaptic connection was determined by a one-to-one 
mapping between the cdf of Ts and that of the target lognormal distribution 
for Vjj (right). (C) Reciprocal and triangle connections are schematically 
shown before (left) and after rewiring (right). In both cases, one of the strong 
connections are chosen randomly and eliminated. 
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Excitatory synapses fail to transmit presynaptic spikes to post- 
synaptic neurons with a certain probability. We denned the failure 
probability as V a /(V a + Vj n ) in terms of the threshold potential 
V a for a synapse with an EPSP amplitude of V,„. Each neuron out- 
puts a spike when its membrane potential reaches threshold Vth, 
then the membrane potential is reset to Vi after the refractory 
period x re f. 

We used the following values of parameters in the present 
numerical simulations. The number of excitatory neurons Ne = 
10,000 and that of inhibitory neurons Ni = 2000. The connec- 
tion probabilities, i.e., the probabilities of finding a non-vanishing 
element in the adjacent matrices, are eg = 0.1 and c\ = 0.5 for 
excitatory and inhibitory synapses, respectively. The decay time 
constants are = 20.0 [ms] and x] n = 10.0 [ms] for excitatory 
and inhibitory membrane decay, and t s = 2.0 [ms] for synap- 
tic decay. We chose average delay of synaptic inputs as x E = 2.0 
[ms] for excitatory synaptic connections, and tA = 1.0 [ms] for 
inhibitory connections. The refractory period is x re f = 1.0 [ms]. 
For membrane potential parameters, spike threshold is V t h = 
—50.0 [mV], reversal potential of leak current is Vi = —70.0 
[mV], and reversal potential of postsynaptic currents are Ve = 
—0.0 [mV] and Vi = —80.0 [mV] for excitatory and inhibitory, 
respectively. Threshold for synaptic transmission failure is V a = 
0.1 [mV]. The weights of synaptic inputs are G EI = 0.017, G IE = 
0.0018, G 11 = 0.0025 for E-to-I, I-to-E, I-to-I connections. For 
E-to-E connections, parameters of lognormal distribution are 
cr = 1.0, and (jl = a 2 + log0.2, and upper limit of EPSP is V max = 
20.0 [mV]. Computer codes are written in C++ and dynamical 
equations are solved by using Euler methods with a time step of 
h = 0.01 [ms]. 

SYNAPTIC CONNECTIONS 

We embedded mutually-independent random memory patterns 
into E-to-E connections as follows. First, we created random 
binary patterns of 0 and 1 {^f },•]_ ! 2 n e accor d m g t0: 



Prob[^ = 1] 



Prob [%? = 0] = 1 



(3) 



where p is the total number of embedded patterns and a is the 
sparseness of these patterns. In this paper, we introduced a con- 
straint as = <zNe to suppress the non-homogeneity across 
different memory patterns. 

The conventional weight matrix of AM network model is 
determined by the local Hebbian rule as • However, 

here we want to create EPSP amplitudes that obey a lognormal 
distribution, while the conventional weight matrix obeys a bino- 
mial distribution. In order to accomplish this, we introduced a 
continuous weight matrix Tu as shown below: 



(4) 



H = l 



where € [0, 1) is a random variable to make the distribution 
of Tjj continuous. The second term in Equation 4 is crucial for 
the genesis of spontaneous activity when the network stores only 
a small number of patterns and hence the first term vanishes 



between most neuron pairs. Then, we determined EPSP Vjj 
between excitatory neurons i and j so that the cumulative fre- 
quency of Vjj may coincide with that of Tjj (Figure IB). Defining 
the set of memory patterns supported by neuron i as = 
{u, \%f = 1 }, we can express Tjj as Tjj = w %f + where 
\Li contains p, = %^ memory patterns. Thus, the cumulative 
function of Tjj is written as: 

5c- 1 

Prob [Tjj <x] = J2 h ifh <?) + (*- x) h (p„ x) 

q = 0 

in terms of h (pj, q) = Pi C q a q (1 — a) p '~ q , where x is the max- 
imum integer below x. On the other hand, the cumulative 
distribution of Vjj is written as: 



Prob [V v <y} = —( 



(logy - u.) 



Z v = i(l + erf[i(logV mflx -n)]). 

where V max is the upper bound for Vjj. Therefore, the map x — >■ y 
is given as follows: 



y : 



exp L + 4la err 1 (2Z,,Prob [x < Tjj] - l) j , 



(5) 



for the pair of neurons i and j. 

If some excitatory neurons send too strong excitatory out- 
puts to other neurons, these neurons can be a potential hazard to 
the stability of network dynamics. To prevent the appearance of 
such hubs, we normalized each synaptic weight over presynaptic 
neurons as: 



y/Zj, Zj = exp ((pj - pa) /pa) , 



(6) 



where the normalization factor Zj is greater or smaller than unity 
if the number of memory patterns encoded by neuron j is larger 
or smaller than the expectation value, respectively. 

As in many other AM network models, the connection matrix 
constructed above is symmetric. In networks of binary or mono- 
tonic analog neurons, symmetric synaptic connections induce a 
"down-hill" motion in the time evolution of network states, mak- 
ing the retrieval of memory patterns possible (Hopfield, 1982, 
1984). The connection matrix shown in Equation 6 generates a 
small number of strong (typically, EPSP >3 mV) reciprocal con- 
nections and triangle motifs. Though the over-representations of 
such motifs have been known in cortical circuits (Song et al., 
2005; Perin et al., 2011), here such motifs form strong exci- 
tatory loops due to the symmetry of the connection matrix. 
Therefore, a small number of neurons belonging to the motifs 
are activated very strongly and tend to burst at very high fre- 
quencies. To avoid this unrealistic situation, we weakened one 
of the edges in all strong loops (Figure 1C). To be precise, we 
searched all "strong" reciprocal excitatory connections that gen- 
erate EPSPs larger than a threshold value V u in both directions. 
Then, we selected one of them randomly and assigned to it 
a new EPSP that is smaller than V\. In addition, we rewired 
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"strong" triangle connections in a similar way such that they 
may not form a circular triangle because neurons in the closed 
loop structure tend to give high-frequency repetitive bursts in 
numerical simulations. In this rewiring, the lower limit of strong 
EPSP Vi is 3.0 [mV] and the upper limit of weak EPSP V u 
is 1.0 [mV]. 

NUMERICAL SIMULATIONS 

We conducted numerical simulations of the model for various 
values of parameters. To initiate spontaneous activity, we applied 
external Poisson spike trains at 10 Hz to all neurons for the initial 
100 ms of each simulation trial and evoked EPSPs with the ampli- 
tude of about 10 mV. Once spontaneous activity is triggered, the 
network continues to show autonomous firing without any exter- 
nal input. This spontaneous firing state may be regarded as the 
resting state of the network. Depending on the values of parame- 
ters, the stability of spontaneous activity is not robust enough and 
the network may exhibit a further spontaneous transition from 
the resting state due to the intrinsic noise generated by reverber- 
ating synaptic input. As explained in the Results, the lifetime of 
spontaneous activity depends significantly on the memory load. 
If the network remained in the resting state for more than 500 ms 
after the initiation of spontaneous activity, we stimulated neurons 
encoding a memory pattern with external input of the duration 
10 ms and the frequency 100 Hz. We may regard this input as a 
cue signal for memory retrieval. 

PATTERN OVERLAPS 

We define a macroscopic order parameter to measure overlaps 
between the network state and memory patterns as: 

i 

where r E = (1/Nb) J2i r > and V = ( l / aN E) %fn> and r > ls 
the firing rate of the z'-th excitatory neuron. In numerical simula- 
tions, the firing rate was calculated from spike data in the interval 
of 800 < t < 1100 [ms]. Spike data during 610 < t < 800 [ms] 
was not used because the system is in a transient state. The pat- 
tern overlap with the retrieved memory pattern maybe called "the 
retrieval rate." 

The retrieval rate is written as K = 1 — r^/r^ (this expression 
is exact due to the restriction on the total number of non- 
vanishing components in each pattern), therefore, K = 1 when 
the system is in the retrieval state and K = 0 in the resting state. 
The retrieval rate K = 0 when (a) all neurons are completely 
silent, (b) when a non-target memory pattern is retrieved or 
(c) when the average firing rate of the retrieval state becomes 
extremely high. The case (a) sometimes occurs when memory 
patterns are not sparse enough (typically, a > 0.15) and the num- 
ber of stored patterns is large. The case (b) occurs regardless of 
the sparseness parameter a when a small number of patterns are 
stored. Note that the target pattern is not stable in this case. In case 
(c), we regarded memory retrieval as unsuccessful since PR neu- 
rons fired at unrealistically high firing rates (typically, >70Hz). 
Therefore, we set K = 0 without calculating the right-hand side of 
Equation (7) although there could be a different treatment for this 



case. Case (c) typically occurs when memory patterns are sparse 
(a < 0.10) and the number of stored patterns is relatively small, 
that is, in the lower left off-diagonal region of the parameter space 
shown in Figures IB and 7B. 

MEAN-FIELD APPROXIMATION 

We analyze the retrieval of a representative memory pattern and 
regard the other patterns as a Gaussian white noise reservoir. 
Then, in this approximation the model consists of three neu- 
ronal populations, i.e., excitatory pattern-retrieval (PR) neurons, 
excitatory background (BG) neurons and inhibitory neurons 
(Figure 2A). In the equations below, indices p and b stand for 
PR and BG, respectively. Below, we explain the outline of the 
mean-field approximation and show the details of the analysis 
in Appendix. This approximation is valid when the total number 
of inputs to each neuron is sufficiently large and the total input 
is approximated to be the sum of a drift term and a fluctuation 
term. In addition, in order for a statistical approach to be valid, 
each memory pattern should contain a sufficiently large number 
of member neurons. 
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FIGURE 2 | Mean-field analysis of the recurrent neural network. (A) The 

approximate treatment employed in the present analysis is schematically 
illustrated. (B) The relation between r E and r s in the balanced state. The 
blue line is the linear regression of simulation results. Each data point (filled 
circle) shows the average firing rates of excitatory and inhibitory neurons in 
both spontaneous state and retrieval state at various values of the 
parameters a and p. 
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We first calculate the average synaptic input from PR neurons 
to PR neurons <G pp > as: 



relationship (Figure 2B). Then, the time evolution of excitatory 
neurons is described by LIF neurons as follows: 



J PPI 



e e *m 



V a + V,j 



-EE 



(8) 



where Fp is a binary probability variable that is defined as 
Prob [F p (X) = l] = X, Prob [F p (X) = 0] = 1 - X. The aver- 
age synaptic inputs <Gbp>, <G p b> and <Gbb> from BG to 
PR, PR to BG, and BG to BG, respectively, are also obtained in 

similar manners. The average squared synaptic weights (g 2 ^ are 
calculated in the same manner. Since <G pp > is larger than the 
average synaptic input averaged over all neurons <Gee>, strong 
recurrent connections among PR neurons contribute significantly 
to sustaining network activity for the recalled memory pattern. 
<Gbp> is typically smaller than <Gee> since each excitatory 
neuron receives almost the same average input. Other connec- 
tions <Gpb> and <G bb > are roughly the same as <Gee> in the 
region of parameter values under consideration. 

If the synaptic weight matrix was a linear superposition of 
memory patterns, G y a Xu, we could separate the con- 

tribution of retrieved pattern and that of non-retrieved patterns 
to the mean field equation and analytically evaluate the storage 
capacity. In our model, however, this separation is quite diffi- 
cult due to a non-linear transformation from Ty = X^ 
to the lognormal synaptic weight matrix G». Nonetheless, we 
can qualitatively understand how the average weight depends 
on the number of stored memory patterns. Noting that the 
synaptic conductance G| £ is approximately proportional to 
the EPSP amplitude which obeys a lognormal distribution, 
we can approximately express the synaptic weight matrix as 
Gf ~ (y/Z) exp(P for infinitely large p, where p = 

a I J pa 2 (l — a 2 ) andZ = exp (Ppfl 2 — (jl)- Thus, we obtain: 



J PPI 



Mi 



-EE 



or < Gpp > I < Gfj >cx exp(l/ v /p), which implies that the rel- 
ative magnitude of reverberating synaptic input among PR neu- 
rons decreases with an increase in the number of stored patterns. 
Although the above approximation is not accurate for finite values 
of p, the asymptotic behavior of <G pp > explains how the storage 
capacity is determined in the present model. 

Then, denoting the average firing rate of each neural pop- 
ulation as Yp, r b , and r,, we obtain the following approximate 
relationship in the balanced state: 



To + k„rE, 



(9) 



where te = ar p + [l-a)r b is the average firing rate of all excitatory 
neurons. Numerical simulations confirmed the validity of this 



l_ 

dt 



, Eq Eq Ne 

^ = -^+E#^E^-*;-^) 



V ') 



j IE IE N, 



y y 



where q = p or b. 

As in the previous model (Teramae et al., 2012), we treated 
excitatory synaptic inputs to each neuron as the sum of Gaussian 
noise generated by weak-dense synaptic inputs and m sparse- 
strong synaptic inputs: 



dt 



+ s E %(t), 



where is a normalized Gaussian random variable with mean 
zero. Then the average and variance of the synaptic conductance 
gf' 1 can be written as: 

{g Ep } = t 5 c E N E (a (Gpp) r p + (l-a) (G bp ) r b ) , 
(g Eb ) = t 5 c E N E (a{G pb ) r p +(l- a) (G b b) ffc) 

= t s c E N E (Gee) m, (10) 
(^) 2 = c E N E (a {G 2 pp) r p + (l-a) {c 2 bp ) r b ) , 

(s^) 2 = c E N E (a (G 2 pb ) rp + a-a) [G 2 hh ) r b ) 



ceNe \Gvv) te 



(11) 



In Equations (10) and (11), we assumed that average inputs to BG 
neurons from PR and BG neurons are approximately the same for 
the two presynaptic neuron categories. Similarly, the fluctuating 
inhibitory synaptic input g 1 ^ is approximated as: 



JE _ IJE\ i m lE dt\¥ ^ T|- £ 



dt 



+ Ai (o ■ 



where (gp) = t s ciN[GiEri, (s 57 ) 2 = C[N[G 2 E r[ and r\ IE is Gaussian 
white noise. Therefore, the membrane potential v q obeys 



ri v l v q - V q 

1 - -L^-^-^-nfevJ-Vj), 



dt 



^ = [l/t E m + (g IE ) + (g Eq )]-\ 

v q = 4[v L /tl + [g Eq )VE + {g IE )Vi\. 
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In the right-hand side of the first equation, we replaced the 
membrane potential with its average value. For this approxima- 
tion to be valid, the variance of the membrane potential of exci- 
tatory neurons should be smaller than Vg — V E and V„ — V[\. 
Otherwise, the fluctuation term may depend on the membrane 
potential and the analytic calculation of firnig rate will be dif- 
ficult. The above equations represent Kramer's equation with a 
linear force. Taking the contributions from strong-sparse synapses 
r 2 , into account, we obtain equations for the time evolution of 
rp and rj, as follows (Risken, 1996; Brunei, 2000; Teramae et al., 
2012): 



dr q 
dt 



(r q -f q (r p j b )). 



r 1 + r 2 . 



■. exp 



(v th -v!) 2 



2o-2 



(12) 
(13) 

(14) 



N s r E 



-erf 



erf 



(v th - vf) 



V2a, 



(Vth - vl - v s ) 



(15) 



where q = p or b, and the variances are given as: 

a 2 Eq = {^)\v q 0 -V E )\ 

aj q = (^) 2 (V^-V I f. 

By solving Equations (12)— ( 15), we derived the various properties 
of the present model reported in the main text. The approxima- 
tion by Kramer's equations is valid when the average membrane 
potential is sufficiently low and the probability distribution of the 
membrane potential over the threshold is negligibly small. The 
approximation is also not so accurate when the average firing 
rate of excitatory neurons is high since the membrane poten- 
tial distribution is biased toward the threshold and spike trains 
are less irregular (in the cases of p = 120 and 125 in Figure 5E). 
We employed a linear approximation (Equation 9) for inhibitory 
neurons since their membrane potentials have a non-negligible 
probability density over the threshold and Kramer's equation is 
not accurate. 

RESULTS 

We conducted numerical simulations of the present model and 
explored its memory retrieval dynamics numerically and ana- 
lytically. In numerical simulations, one of the stored patterns 
(t was excited by external input. We found that a heavy-tailed 
weight distribution of Hebbian synapses creates memory- specific 
subthreshold membrane potential fluctuations to regulate the 
stochastic dynamics of memory retrieval. In the rest of the paper, 
we call excitatory neurons encoding the retrieved pattern as PR 



(pattern retrieval) neurons, and those not encoding the pattern 
as BG neurons. For instance, when pattern 1 is retrieved (u, = 1), 
neuron i is a PR neuron if ^! = 1 or a BG neuron if ^! =0. 

Figure 3A shows an example of successful memory retrieval by 
the network model. Without external input, the model is able to 
sustain spontaneous activity in which both PR neurons and BG 
neurons fire at low firing rates. If PR neurons (encoding mem- 
ory pattern 1) are innervated by a brief cue signal (Materials and 
Methods), the model changes its dynamic behavior from sponta- 
neous activity to a retrieval state, in which the average firing rates 
of the corresponding PR neurons are much higher than those of 
BG neurons. In fact, the firing rates of BG neurons are lower than 
the average firing rate of all excitatory neurons. As explained later, 
the network model can retrieve a memory pattern successfully if 
the number of stored patterns is within a certain range between 
upper and lower critical values. When this number is close to 
the upper bound (i.e., the storage capacity), excitatory neurons 
show sparse irregular firing in the retrieval state (Figures 3B,C). 
Inhibitory neurons also exhibit high coefficients of variations 
(CVs), although they fire near-synchronously at relatively high 
firing rates (Figure 3D). 

When the number of stored patterns is smaller than a lower 
critical value, the model does not have a stable spontaneous activ- 
ity. Due to the intrinsic noise generated by recurrent synaptic 
input, the network state eventually evolves from the resting state 
into one of the embedded patterns (Figure 3E). Therefore, the 
resting state is only quasi-stable. Furthermore, the memory pat- 
terns evoked in the final network state depend neither on the 
initial state nor cue signal. This is presumably because the state 
transition almost always targets such memory patterns that are 
separated by lower potential barriers from the resting state than 
other patterns. Therefore, the other memory patterns are difficult 
to recall and the model is unable to perform AM. If p is close 
to, but larger than the lower critical value, the network model 
remains in the initial resting state for more than hundreds of mil- 
liseconds, much longer than the membrane time constant. This 
long lifetime of the resting state indicates that the state transition 
is caused by stochastic fluctuations in network activity. Even if 
the spontaneous activity state no longer exists, the network could 
show slow dynamics due to the flatness of the potential func- 
tion around the ghost of spontaneous activity state. In this case, 
the lifetime of the (ghost) spontaneous state would not be dis- 
tributed exponentially. In practice, it is difficult to find the exact 
critical value by numerical simulations or by the present mean- 
field approximation, which does not explicitly take into account 
the influences of non-retrieved memory patterns. 

Regarding the trial average <K> of the parameter K, which 
was evaluated over 15 trials on each value of parameters, as 
the success rate of retrieval, we display the areas of the two- 
dimensional parameter space spanned by p and a in which <K> 
is greater than 50% (Figure 3F). Given the value of a, the net- 
work model can retrieve a memory pattern successfully only if 
the number of stored patterns is between upper and lower criti- 
cal values, where the upper bound represents the critical storage 
capacity (Hopfield, 1982; Amit et al, 1985). Interestingly, the 
retrieval states become unstable also in the regime of low memory 
load. This property is characteristic to the models that replicate 
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FIGURE 3 | Simulation of memory retrieval in the associative memory 
model. Here, a = 0.12 and p = 140. (A) Raster plot and mean firing rates are 
shown. In the raster plot, red dots represent spikes of inhibitory neurons and 
black ones are spikes of excitatory neurons. PR neurons were displayed at 
the bottom. Blue vertical bars shows a period in which PR neurons were 
stimulated selectively by external inputs. The firing rates were calculated 
with a time bin of 10 ms. Gray lines are the firing rates of neurons encoding 
(p-1) non-retrieved memory patterns. (B,C) The distributions of firing rates 
and CVs for PR neurons and inhibitory neurons. (D) The power spectrum of 



inhibitory neural activity. We used the firing rates that were calculated with 
time bins of 1 ms. The results shown in (B), (C), and (D) were calculated from 
the simulated spike data in the interval 650 ms <f<1600 ms shown in (A). 
The average firing rates of BG neurons are too low (0.9 Hz) to calculate the 
distributions of firing rates and CVs. (E) Simulation results are shown at 
a = 0.12 and p = 110. Note that the value of p is smaller than the one used in 
(A), and one of memory patterns was activated without a cue signal around 
300 ms. (F) The average success rate <K> was numerically obtained for 
various values of p and a. 



spontaneous activity (Curti et al., 2004; Latham and Nirenberg, 
2004; Roudi and Latham, 2007), and is not seen in other t\M 
networks that do not have such a global state. An intriguing 
property of our model is that it generates spontaneous activity 
internally by reverberating synaptic input without external noise. 
We will demonstrate the non-trivial effect of internal noise later 
in detail. 

ROLE OF THE SUBTHRESHOLD MEMBRANE POTENTIALS IN MEMORY 
RECALL 

Long-tailed distributions of EPSPs were previously shown 
to achieve stochastic resonance between postsynaptic firing 



and presynaptic spikes at sparse strong recurrent synapses in 
asynchronous irregular states of cortical networks (Teramae et al., 
2012). In this stochastic resonance, the average membrane poten- 
tials of individual neurons are maintained at a subthreshold level 
optimal for spike transmission by very strong synapses. Below, we 
demonstrate how the present network controls the subthreshold 
membrane potential dynamics of individual neurons for efficient 
memory recall. 

In the parameter regime of the bistable network dynamics, the 
average membrane potentials of PR and BG neurons split into two 
distinct levels in the retrieval state: the average membrane poten- 
tials are more depolarized in PR neurons than in BG neurons 
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FIGURE 4 | Role of the subthreshold membrane potential distributions 
in memory retrieval. (A) Time evolution of the average membrane 
potential is shown for a = 0.12 and p= 140. Blue vertical bar shows the 
duration of stimulus presentation as in Figure 3A. (B) Relationship 
between the average membrane potential and the number of stored 
patterns at a =0.12. Gray squares show the average membrane potential 
of all excitatory neurons in spontaneous activity state, whereas blue and 
green squares are the averages of PR neurons and BG neurons in 
retrieval states. Error bars show the standard deviations of the trial 



variability. (C) Distributions of the membrane potential were calculated for 
spontaneous and retrieval states. In calculating the distributions for BG 
neurons, N E a neurons were randomly chosen for the normalization. Peaks 
at v = —70 mV, which arose from the refractory period, were removed 
from the distributions. (D) Relationships between the spiking probability 
and average membrane potential. We calculated the spiking probability of 
a postsynaptic neuron from 1 to 5 ms after the firing of the presynaptic 
neurons sending strong inputs (EPSP >6mV) to the post-synaptic 
neuron. 



(Figure 4A). The differences in the average membrane potentials 
between the two categories of neurons become smaller as the 
memory load becomes heavier (Figure 4B). In the spontaneous 
firing state, distinctions between PR neurons and BG neurons are 
merely formal and they have identical Gaussian-like membrane 
potential distributions with the same means and variances. In 
the retrieval state, the membrane potentials of PR and BG neu- 
rons also obey Gaussian-like distributions with approximately the 
same variances. However, their means are clearly different for the 
two classes of neurons (Figure 4C). This splitting of the EPSP- 
amplitude distributions has significant implications in memory 
retrieval. The amplitude of EPSPs is typically larger than 6 mV 
at very strong synapses. This means that 10% of neurons can 
fire in the spontaneous firing state when a presynaptic spike 
arrives at a very strong synapse. By contrast, in the retrieval 
state the membrane potential distributions are shifted toward 
a more depolarized level in PR neurons (v, > — 56 mV), hence 
their probability of firing in response to such a presynaptic input 
increases by more than 30%. On the contrary, BG neurons rarely 
have membrane potentials higher than — 56 mV, so they show 
only a very small firing probability for input at a strong synapse 
(Figure 4D). 



COMPARISON WITH THE MEAN-FIELD APPROXIMATION 

To clarify the dynamical properties of the present model, we 
conducted an analytical study using the mean-field approxima- 
tion (Materials and Methods) and examined the validity of the 
analytical results by numerical simulations. In the mean-field 
approximation, we deal with the retrieval of a representative 
memory pattern by regarding other patterns as a noise reservoir. 
In this approximation, the model consisted of the three neuronal 
populations, i.e., excitatory PR neurons, excitatory BG neurons, 
and inhibitory neurons (Figure 2 A). Below, the indices p and b 
refer to any quantity related to PR or BG neurons, respectively. 

To investigate the dynamical phases of the model, we derived 
the nullclines of the average firing rates of PR and BG neurons 
by setting as r* = ry = 0 in Equations (12)— (15). In doing so, 
we assumed that inhibitory neurons are enslaved by excitatory 
neurons. In deriving the nullclines, we used parameters rj„ and 
k„ to fit the balanced firing rates obtained in numerical simula- 
tions (Figure 2B) and derived the averages and variances \Gqd) 

and (<J, 4 = P> ^> E) from numerically generated Hebbian 

connection matrices G«. We evaluated the contributions of the 
long tail of synaptic connections to the mean-field equations by 
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taking into account the contributions of N s strong-sparse connec- 
tions and those of weak-dense synapses separately. This number 
remains to be a free parameter that should be fixed at a reason- 
able value. We found that the choice of N s = 3 and the EPSP 
amplitude of V s = 6.5 mV for these synapses yield an excellent 
agreement between the analytical and numerical results. 

When the number of stored patterns is in an adequate range, 
the nullclines of r p and r\, yield three intersection points cor- 
responding to the fixed points of the network dynamics in the 
two-dimensional space spanned by the parameters. According to 
a linear stability analysis, the central fixed point is unstable and 
the two fixed points on both sides are stable. Therefore, a stable 
fixed point in the regime of low or high r p corresponds to spon- 
taneous activity and the retrieval state, respectively (Figure 5 A). 
As the number of stored patterns is increased, the retrieval state 
disappears at a critical value through a saddle-node bifurcation 
(Figure 5B). As a result, the system has only one stable state 
corresponding to spontaneous activity. 



We can also determine the storage capacity by numerical sim- 
ulations. As more patterns are stored the value of <K> gradually 
decreases until it finally vanishes (Figure 5C). The value of <K> 
is close to 0.5 in the vicinity of the critical point. Therefore, we 
may define the storage capacity as the number of patterns for 
which <K> is 0.5. The values of the storage capacity thus eval- 
uated numerically and analytically are shown in Figure 5D as a 
function of the sparseness parameter. In Figure 5E, we also calcu- 
lated the average firing rates of spontaneous activity rg and those 
of PR and BG neurons in the retrieval state for a = 0.12. In all 
cases, the numerical and analytical treatments show a reasonably 
good agreement. 

NETWORK DYNAMICS AT LOW MEMORY LOAD 

As mentioned previously, this model also has a lower bound for 
the memory load in addition to an upper bound or the storage 
capacity. This instability reflects the fact that the model gener- 
ates internal noise for memory recall by weak-dense synapses. 
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FIGURE 5 | The mean-field analysis of the equilibrium states. (A,B) The 

nullclines of r p (red) and r b (blue) are shown for the values of p that are 
well below or greater than the critical storage capacity. Filled and empty 
circles show stable and unstable fixed points, respectively. (C) Relationship 
between <K> and the number of stored patterns is plotted. Error bars 
represent SD. We calculated the average value of of over different 
patterns \i in a given network and then averaged the value over different 
realizations of the network (i.e., the connection matrix). Therefore, 



M = E k k m Ti m " Kf'AWu™*) and the standard deviation 

Cf= ^/k max yjj™* (l/Hmax EH™* 4 k) - (K)f ■ where index k runs over 
the different realizations and k max = \imax = 10. (D) The storage capacity 
was evaluated numerically and analytically as a function of the sparseness 
parameter. (E) The average firing rates of all excitatory neurons were 
calculated in various states at a = 0.12. "SA" labels spontaneous 
activity. 
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Below, we investigate this instability by means of the mean-field 
approximation. 

Because fluctuations in r/, are much smaller than those in r p , 
we may replace rj with its fixed point r? in the mean-field analysis. 
Then, the dynamics of the model is approximately described with 
a single parameter r p as: 



dr 



_P_ 
dt 



1 

~p{ r P 



(r p - f p (rp, r* b )) 



dU(r p ) 



dr 0 



(16) 



We calculated this potential U (rp) numerically. When the num- 
ber of embedded patterns is in the range that enables memory 
recall, the potential 17 (r p ) has two local minima corresponding 
to spontaneous activity and the retrieval state at r p = s m ;„ and 
r m i„, respectively (Figure 6A). As the number of stored patterns is 
decreased, the latter local minimum becomes deeper (Figure 6A, 
inset) and simultaneously the potential barrier AE separating 
the two local minima becomes lower (Figure 6B). Note that we 
measure AE always from the bottom of the potential, t7(s m ;„), 
in the present argument. The result implies that spontaneous 
activity cannot remain stable when the number of stored patterns 
becomes smaller. In this case, fluctuations in spontaneous activ- 
ity make it easier for the network state to jump over the potential 
barrier into the retrieval state. On the contrary, when the num- 
ber of stored patterns exceeds a critical value, the local minimum 
corresponding to the retrieval state disappears and the network 
ceases to function as AM (Figure 6C). Figure 6D displays the 
relationship between AE calculated by the mean-field approxi- 
mation and the mean lifetime of spontaneous activity obtained 



by numerical simulations. The mean lifetime depends exponen- 
tially on AE, suggesting that the transition to the retrieval state is 
due to stochastic fluctuations around the spontaneous firing state. 

While this model accomplishes a relatively high storage capac- 
ity with spiking neurons, it does not show a stable and robust 
retrieval of an arbitrary memorized activity when the number 
of stored memory patterns is small. While long-tailed EPSP dis- 
tributions are known in neocortical and hippocampal circuits 
(Song et al, 2005; Lefort et al, 2009; Ikegaya et al, 2013), the 
instability of the model seems to be unrealistic in biological 
neural networks. Therefore, we explored a way to modify the 
model for circumventing this difficulty. We point out that the 
unrealistic property arises from the theoretical hypothesis that 
the weights of all excitatory synapses are determined solely by 
stored memory patterns. Therefore, we created a partly random- 
ized connection matrix and studied the retrieval dynamics of 
a neural network with such synaptic connections by numerical 
simulations. 

We define partly randomized EPSP amplitudes V| r of E-to-E 
synapses as follows: 



n 



v y - ~ logN([i,a) 

Va 



ifr|, ; < c r 
otherwise 



(17) 



where Vy is the EPSP amplitude constructed previously from the 
long-tailed EPSP distribution, c r is the fraction of the randomized 
synapses in all E-to-E synapses (0 < c r < 1), and r|y is a random 
variable that takes a value between 0 and 1 (Figure 7A). Here, c r 
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FIGURE 6 | The landscape of the potential fucntion U[r p ). (A) The 

potential function calculated at a = 0.12 and p= 135 has two local minima 
at r p = s m / n and r mjn corresponding to the spontaneous activity and 
retrieval state in Figure 5A, respectively. Inset displays a similar potential 
function for a smaller number of memory patterns (a = 0.12 and p= 120). 
(B) Relationships between the number of stored patterns and potential 
barrier AE that separates the spontaneous acivity state and the retrieval 



state for different values of sparseness. Here, p= 135. (C) The potential 
function calculated at a = 0.12 and p= 150 has only a single minimum 
corresponding to the stable fixed point in Figure 5B. (D) Relationship 
between the analytically calculated potential barrier and the mean life time 
obtained in simulations. Here, a =0.1. The mean life time is the time the 
network took for escaping from initial resting states. Error bars 
show SD. 
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FIGURE 7 | Associative memory networks with Hebbian and random 
synaptic components. (A) The distribution of the random components in 
synaptic weights for c r = 0.4. The portion of synaptic connections used for 
storing memory patterns are colored as in Figure IB, while the gray area 
shows the portion of synaptic connections not employed for the memory 
storage. (B) Phase diagram of the success rate <K> for various values of 
c, and the number of embbed memory patterns p. (C) Raster plot and the 
time evolution of the averaged firing rates. Parameter values were a = 0.10, 
c, = 0.2, and p = 100. Both panels were drawn in the same way as in 
Figure 3A. 



cannot be too large to maintain reasonably large memory com- 
ponents in the connection matrix. The conductance constants 
were set as G EI = 0.016, G IE = 0.0018, and G 11 = 0.002, and the 
other parameters took the same values as in the previous simu- 
lations. If we increase the noise fraction c r for a small number 
of memory patterns, the retrieval state is also stabilized with a 
small memory load (Figure 7B). Spiking activity of the modified 
model during memory retrieval is similar to that of the previ- 
ous model (Figure 7C). Thus, the inclusion of noisy components 
in synaptic connections secures the stability of retrieval dynam- 
ics with Hebbian synapses obeying long-tailed EPSP amplitude 
distributions. 



DISCUSSION 

We have presented an AM network model of spiking neurons 
and explored its dynamical properties during memory retrieval. 
According to results of recent electrophysiological studies (Song 
et al, 2005; Lefort et al, 2009; Ikegaya et al, 2013), we con- 
structed Hebbian synaptic connetions that obey a long-tailed or a 
lognormal distribution of synaptic weights. Our model possesses 
a spontaneous firing state with low frequency neuronal firing 
(~1 Hz) and multiple persistent states with high frequency neu- 
ronal firing (~30 Hz), in which one of the embedded memory 
patterns is recalled. 

IMPLICATIONS OF LONG-TAILED WEIGHT DISTRIBUTIONS IN 
STOCHASTIC NETWORK DYNAMICS 

Several papers have already studied the implications of long-tailed 
distributions of EPSP amplitudes in the dynamics of local corti- 
cal networks. Such EPSP-amplitude distributions account for a 
long-tailed distribution of firing rates in spontaneous firing of 
in vivo cortical neurons (Koulakov et al., 2009). The distribution 
contains many weak synapses and a small number of extremely 
strong synapses, and we recently showed that this coexistence 
of dense-weak and sparse-strong synapses enables recurrent cor- 
tical networks to generate an autonomous activity with highly 
irregular sparse firing of single neurons (Teramae et al., 2012). 
This persistent activity replicated various statistical properties of 
spontaneous cortical activity observed in experiment. We demon- 
strated that long-tailed weight distributions generate intrinsic 
noise optimal for spike communications by individual neurons. 
In short, the stochastic resonance effects generated by massively 
many weak synapses facilitate the faithful transmission of spike 
information received at strong synapses. Ikegaya et al. (2013) 
recorded long-tailed weight distributions also in the hippocam- 
pus and constructed a recurrent network model based on their 
experimental findings. This model involved NMDA receptors at 
recurrent synapses, which were crucial for stabilizing sponta- 
neous activity in the model. The present model, however, does not 
involve NMDA receptors as the slow synaptic current is unneces- 
sary for the genesis of spontaneous activity (Teramae et al., 2012). 
Finally, we note that long-tailed distributions of synaptic weights 
are well consistent with STDP (Gilson and Fukai, 201 1). 

COMPARISON WITH PREVIOUS SPIKING MODELS OF ASSOCIATIVE 
MEMORY 

Several models have been proposed to study the dynamics of 
memory retrieval in AM networks of integrate-and-fire type 
neurons (Sommer and Wennekers, 2001; Latham and Nirenberg, 
2004; Curti et al, 2004; Aviel et al, 2005; Roudi and Latham, 
2007). These models typically store sparsely coded memory pat- 
terns, for which AM networks in general exhibit a good per- 
formance. In comparison with AM networks of binary neurons, 
however, the networks of spiking neurons can generally store a 
relatively small number of memory patterns (Curti et al., 2004; 
Latham and Nirenberg, 2004; Roudi and Latham, 2007). Though, 
to our knowledge, no rigorous proof has been known, this low 
performance presumably reflects the leaky property of synaptic 
integration by spiking neurons. This is an interesting difference 
between the two classes of neural network models. It is also 
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often difficult to maintain the firing rates of PR and BG firing 
neurons within a physiologically reasonable range, while keeping 
irregular spiking as typically observed in in vivo cortical neurons. 

It is worthwhile to compare the storage capacity between dif- 
ferent AM models of spiking neurons. A regorous comparision 
seems to be difficult since the retrieval performance of such 
models often shows a strongly non-linear dependence on param- 
eter values. However, the storage capacity in general increases 
with the sparseness of memory patterns and decreases with the 
sparseness of synaptic connectivity. Therefore, the parameter 
ct= (pa | In a | /ceNe) may be used for measuring the storage 
capacity (or the information capacity) of AM models with sim- 
ilar sparseness and network size. The value of this parameter is 
0.036 (N E = 10,000, c E = 0.1, a = 0.12, p = 140) in our model, 
0.0058 (N E = 8000, c E = 0.25, a = 0.1, p = 50) in Latham and 
Nirenberg (2004), and 0.0056 (N E = 8000, c E = 0.2, a = 0.05, 
p — 60) in Curti et al. (2004). Therefore, out model has a larger 
storage capacity than the other models of spiking neurons in 
similar ranges of parameter values. 

This model shows spontaneous firing of neurons, or BG activ- 
ity, without any external input. As shown in our previous model 
of spontaneous cortical activity (Teramae et al., 2012), the coex- 
istence of sparse-strong and dense-weak synaptic connections 
on each neuron enables individual neurons to faithfully trans- 
mit presynaptic spikes at sparse-strong synapses by means of 
stochastic resonance. In other words, the model utilized dense- 
weak recurrent synapses to generate internal noise optimal for 
this resonance effect. Input to dense-weak synapses maintains 
a depolarized subthreshold membrane potential and input to 
sparse-strong synapses evoke spikes with a certain probability. 
Thus, the broad range of synaptic strength enables the individual 
neurons to maintain spontaneous firing at low firing rates. 

Our study has further shown that internal noise can selec- 
tively sustain the persistent firing of neuronal subgroups when 
synaptic connections are determined by the Hebbian rule. Due 
to relatively strong recurrent connections between PR neurons, 
they receive intense noise via weak-dense synapses when they are 
activated. This noise strongly depolarizes their membrane poten- 
tials and consequently these neurons fire with high probabilities 
in response to a spike input at sparse-strong synapses. By con- 
trast, BG neurons have lower subthreshold membrane potentials 



and lower probabilities of spike responses to strong inputs. This 
probabilistic spiking achieved by the cooperation of SSWD synap- 
tic inputs may provide a novel mechanism of probabilistic neural 
computations. 

It was previously pointed out that spontaneous firing states 
turn to be unstable in AM network models of spiking neurons 
when the number of stored memory patterns is small (Curti 
et al., 2004). It is also known that when synaptic connections are 
described as the sum of Hebbian components and random com- 
ponents, the Hebbian components should be sufficiently small to 
maintain stable spontaneous firing states (Latham and Nirenberg, 
2004; Roudi and Latham, 2007). The present model also exhibits 
the instability of spontaneous firing states when the memory load 
is too low. According to the previous results, we demonstrated 
that the addition of random components to lognormally dis- 
tributed Hebbian components in the connection matrix removes 
the instability from the model. 

The genesis of highly irregular asynchronous states is not 
trivially easy in modeling persistent activity of recurrent neu- 
ral networks. Such activity is typically seen in working mem- 
ory, and various solutions to this problem have been proposed, 
which include the uses of slow reverberating synaptic current 
(Wang, 1999), balanced excitatory and inhibitory synaptic input 
(Renart et al, 2007; Roudi and Latham, 2007; Mongillo et al, 
2012), short-term synaptic depression (Barbieri and Brunei, 2007; 
Mongillo et al., 2012), and modular network structure (Lundqvist 
et al., 2010). In this study, we have shown that long-tailed distri- 
butions of EPSP amplitudes generate a highly irregular persistent 
activity in the memory retrieval, as was the case in spontaneous 
activity (Teramae et al, 2012). 

In summary, we have developed an AM model of spiking neu- 
rons by taking long-tailed distributions of synaptic weights into 
account. Our network model exhibited excellent performance in 
the memory retrieval, suggesting an active role of internal noise 
in memory processing by the brain. 
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APPENDIX 

DERIVATION OF FLUCTUATIONS IN SYNAPTIC INPUTS s Ep , s Eb 

We define <X> as the averaging over both neural population and time bins. Amplitude fluctuations in excitatory synaptic inputs 
on PR neurons are represented as: 



= ( E C 'J G 7 E h(t ~ S J ~ d ^ ~ CeNe ( a < G PP >r P + (l-a)< G bp > r b ))- 



■■jCff ^2 8 ( f - s i ~ d v) - a < G pp > r p 



+ 



E 



q.jG? ^ 8(t - Sj - dij) - a < Gbp > 



rb 



Thus, they can be separated into two components corresponding to input from PR neurons and that from BG neurons. For the former 
input, 



CijCff ^ 8(f - s,- - dij) - a <G pp > r p 



4 (of- < Gpp >) \J2 W ~ s i ~ d 'i^ 



p 

E 

j 



+ 



E 



4 < Gpp > 2 E 8 ^ _ S J ~ ^ ~ r P 



N E a [cg(l - c E ) < G 2 pp > r ? +c| (< G 2 pp > - < G pp > 2 ) r p + c\ < G pp > 2 rp] 



(Al) 



(A2) 

c E N E a < Gpp > r p . (A3) 
In the calculation of Equation (A3), we assume that the input is in the regime of low firing rates. Therefore, {^2 Sj 8 (f — Sj — d,j 



0 or 1, and 



E>-9 ^(^Ht-sj-di^. 



(A4) 



Furthermore, we may assume that synaptic input is represented by Poisson spike trains. Under this assumption, the fluctuation of 
firing rates are written as, 



I^W-Sj-fy-rA 



Input from BG neurons can be similarly calculated and the variance is evaluated as: 

(s^) 2 = c E N E (a {G 2 pp) r p + (1 - a) (g^) r b ) , 

and (s Eb ) 2 is evaluated in a similar manner. 
DERIVATION OF THE POPULATION RATE DYNAMICS 

In Gaussian approximation of synaptic inputs, the membrane potential v q is written as follows: 



(A5) 



(A6) 



d± 
dt 



\ {A - O - nf g iyl - v E ) - nf (v* - v,) . 



4=[l/r E m + ^ E ) + (^)]-\ 

v! = xl[V L /r E m + { g ^)V E + {g IE )V I ]. 

In the calculations below, we drop index i from equations as they are essentially the same for all excitatory neurons. 
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Let vj = y q , —r\f q (V<j — Ve) = o£, —\ q (v| — Vj) = xl, then the above equations are rewritten as: 

where (ct^)' = (s E if (V q 0 - V E f , (otf = (/«) 2 (V q 0 - V;) 2 . Defining = -| + x\ + x£, we obtain 

^9 _ 



i i 

^7 + - 



z-^- + JK) 2 + (^) 2 ? 9 a). 



Using Kramer's equation, we can derive the stationary probabilistic distribution W(y, z) for the above two variables as: 



W(y, z) 



27TO-2 



exp 



t e t5 2 _ r_ 



1 1 

4 + ^ s 



1 xUs 



Then, we can calculate the output firing-rate of each neuron from the probabilistic current as: 



f 

Jo 



zW(V th - V„, z)dz : 



1 



2jt*/x e x s 



exp 



Mr 



(A7) 



An additional cause of the rate increases in postsynaptic neurons is the pool of extremely strong EPSPs, which contribute to the firing 
rate of each neuron by: 



N s r E 



erf 



v/2c 



(V th - V q ) 



erf 



\/2a. 



{v th - vl - v,) 



(A8) 



where JV S is the number of extremely strong synapses and V s is the average EPSP size of these strong inputs. Thus, the overall firing 

« = r l + r \ 



rate is given as r q = r\ + r ( 2 
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